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O : ABSTRACT 

S ' Non-relativistic three-dimensional magnetohydrodynamical (MHD) simulations 

O ■ of Poynting flux dominated (PFD) jets are presented. Our study focuses on the prop- 

^ ■ agation of strongly magnetized hypersonic, but sub-Alfvenic (C^ ^ V-^^ < V^) flow 

I ■ and on the subsequent development of a current-driven (CD) kink instability. This 

;_( ■ instability may be responsible for the "wiggled" structures seen in sub-parsec scale 

^ ■ (VLBI) jets. In the present paper, we investigate the nonlinear behavior of PFD jets 

in a variety of external ambient magnetized gas distributions, including those with 

1^. density, pressure, and temperature gradients. Our numerical results show that the jets 

'j_j ■ can develop CD distortions in the trans-Alfvenic flow case, even when the flow itself 



is still strongly magnetically dominated. An internal non-axisymmetric body mode 
grows on time scales of order of the Alfven crossing time and distorts the structure 
and magnetic configuration of the jet. The kink (m = 1) mode of the CD instability, 
driven by the radial component of the Lorentz force, grows faster than other higher 
order modes (m > 1). In the jet frame the mode grows locally and expands radially 
at each axial position where the jet is unstable: the instability, therefore, does not 
propagate as a wave along the jet length. CD instabilities have a number of features 
that make them an attractive explanation for the helical jet structure observed in AGN 
and pulsars: 1) because the magnetic field remains strong, CD instabilities do not de- 
velop into full MHD turbulence; 2) the helical structures saturate and advect with the 
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bulk flow; 3) they distort the body of the jet, not merely its interface with the ambient 
medium; 4) local plasma flow, then, follows a helical path along the kinked magnetic 
field backbone. A naturally-occurring, external helically magnetized wind, which is 
(quasi-) axially current-free, surrounds the well-coUimated current-carrying jet and re- 
duces velocity shear between the jet and external medium. This stabilizes the growth 
of MHD Kelvin-Helmholtz surface modes in the inner jet flow. 

Subject headings: Instabilities — galaxies: active — galaxies: jets — methods: numerical- 
magnetohydrodynamics (MHD) 



1. INTRODUCTION 

Magnetohydrodynamic (MHD) acceleration mechanisms often are invoked as a model for 
the launching and initial acceleration and coUimation of winds and jet outflows from Young Stel- 
lar Objects (YSOs), X-ray binaries (XRBs), Active Galactic Nuclei (AGNs), Microquasars, and 
Quasars (QSOs) (see, e.g., Meier, Koide, & Uchida 2001, and references therein). There has been 
a growing recognition in recent years, however, that the influence of strong magnetic fields within 
the jet may extend well beyond the central engine into the region where the jet freely propagates, 
influenced only by internal and ambient hydro- and magnetohydrodynamic forces. This is partic- 
ularly evident in observations of jets in AGNs, QSOs, winds from pulsars, and y-ray burst sources 
(GRBs) (e.g., Perley, Bridle, & Willis 1984; Conway & Murphy 1993; Hester et al. 2002; Coburn 
&Boggs 2003). 

Strongly magnetized jets, particularly those with a strong toroidal field encircling the coUi- 
mated flow, are often referred to as "current-carrying" or "Poynting flux dominated" (PFD) jets. A 
large current flowing parallel to the jet flow is responsible for generating a strong, tightly-wound 
helical magnetic field. Continued rotation of the entire magnetized plasma about the jet axis also 
plays an important role in jet dynamics. Rotation of the helical field drives a "barber pole"-like, 
torsional Alfven wave (TAW) forward in the direction of the jet flow, carrying electromagnetic 
energy and further accelerating the plasma. In a PFD jet, the Poynting flux energy carried by this 
TAW can greatly exceed the kinetic energy flux in the hydrodynamic (HD) part of the flow. 

It is well-known that a cylindrical plasma column with a helical magnetic configuration is sub- 
ject to MHD instabilities. These are usually divided into pressure-driven (PD), Kelvin-Helmholtz 
(KH), and current-driven (CD) instabilities (see e.g., Kadomtsev 1966; Bateman 1980; Freidberg 
1982, for details). PD instabilities are related to the interplay between gas pressure and curvature 
of magnetic field lines (Kersale, Longaretti, & Pelletier 2000; Longaretti 2003), and have not been 
considered to be very important for a supersonic jet. KH instabilities occur because of the pres- 



ence of velocity gradients in the flow (see e.g.. Landau & Lifshitz 1959; Chandrasekhar 1961, 
for details). They may play an important role at the shearing boundary between the flowing jet 
and external medium, particularly in kinetic flux dominated (KFD) jets where the hydrodynamical 
forces dominate over magnetic. PFD jets, on the other hand, should be especially susceptible to CD 
instabilities because of the presence of the strong axial electric current. The investigation into the 
destructive influence of CD instabilities on jet flow, therefore, recently has become an important 
avenue of research in the study of astrophysical jets. 

The purpose of the present paper is an in-depth, 3-D MHD, non-relativistic numerical inves- 
tigation of the nonlinear development of CD instabilities in PFD jets — particularly the CD kink 
(m = 1) mode. Our jets are endowed with most of the properties now thought to characterize 
those in AGN, QSOs, pulsars, and GRBs: strongly magnetized, hyper-sonic (but sub-Alfvenic; 
C^ ^ V-^ < V^) flow, driven by a Foynting flux dominated, torsional Alfven wave that carries most 
of the energy. In a previous paper (Nakamura, Uchida, & Hirose 2001), the basic behavior of 
TAWs and PFD jets was studied under simplified atmospheric conditions. We now assume more 
realistic atmospheric situations, including density, pressure, magnetic field, and temperature gra- 
dients in the ambient medium. While these jets are strictly ''propagating" (we do not consider the 
central engine itself), it is important to realize that the inclusion of toroidal magnetic fields and 
rotation can create a local acceleration and coUimation process within the jet that can increase the 
speed of the jet downstream up to and beyond the Alfven speed. Unlike previous studies of prop- 
agating jets, therefore, the downstream properties of the jet will be determined not only by the flux 
of kinetic energy injected at the jet throat (i.e., the initial jet speed andMach number), but also by 
the amount of Poynting flux injected in the TAW itself. 

In section 2, we give a review of observational and theoretical work on the subject, relating 
previous results to the current work. Section 3 outlines our basic numerical methods and model. 
Section 4 gives a comprehensive report and discussion on our results. Our conclusions are sum- 
marized in section 5. 



2. MAGNETOHYDRODYNAMIC vs. HYDRODYNAMIC JETS 

2.1. Observational Evidence for Strongly Magnetized Jets 

There are several observational results which indicate the presence of strong magnetic fields 
in astrophysical jets. The first class of observations concerns the structure of the magnetic field. 
In the case of AGN jets, this is determined from radio polarization of the synchrotron emission, 
especially the Rotation Measure (RM) and the projected magnetic field vectors {e.g., Perley, Bridle, 
& Willis 1984; Owen, Hardee, & ComweU 1989; Perlman et al. 1999; Feretti et al. 1999; Eilek 



& Owen 2002; Krause & Lohr 2004). Asada et al. (2002) showed that the RM distribution for 
the quasar 3C 273 jet on parsec scales has a systematic gradient across the jet, and the projected 
magnetic field vector is systematically tilted from the direction of the central axis of jet. They 
conclude that the sign reversal of the RM across jet indicates the presence of a toroidal (azimuthal) 
component of the field inside jet. In the case of y-ray burst (GRB) jets, Coburn & Boggs (2003) 
observed linear polarization in the prompt y-ray emission from GRB021206 and found a value 
of 80 ±20%. This is the theoretical maximum possible polarization for a magnetized plasma, 
indicating that the field is so strong that the entire y-ray-emitting region is organized around the 
field structure. The flow clearly is magnetically dominated, and the authors further suggest that the 
GRB explosion itself is powered by magnetic fields. 

A second class of observations concerns the morphological structure of the jets themselves. 
High resolution VLBI observations show that many AGN and quasar jets display wiggles or kinks 
on sub-parsec to parsec scales (e.g., Krichbaum et al. 1990, 1992; Hummel et al. 1992; Conway 
& Murphy 1993; Roos, Kaastra, & Hummel 1993; Jones et al. 1996; Krichbaum et al. 1998; 
Mantovani et al. 1999; Murphy et al. 1999; Hutchison, Cawthorne, & Gabusda 2001; Lobanov 
& Zensus 2001; Stirling et al. 2003). Such a helical distortion might be caused either by MHD 
instabilities, or by precession of the jet ejection axis due to the existence of a binary Black Hole 
(BBH) (Begelman, Blandford, & Rees 1980), or by an encounter with another galactic core. On 
the other hand, the precession and galaxy interaction models tend to operate on rather long (10^ 
yr) time scales, while the kinks appear to occur on significantly shorter (<^ 10^ yr) time scales 
{e.g.. Potash & Wardle 1980). Based on both theoretical considerations and observational results, 
therefore, MHD instabilities appear to be the most plausible model for the observed wiggled or 
kinked structures. 

The third indicator of possible dynamical importance of magnetic fields is the presence of 
thermal overpressures in jets. There are several such cases in which the jet thermal pressure sig- 
nificantly exceeds the surrounding X-ray gas pressure (see e.g.. Potash & Wardle 1980; Owen, 
Hardee, & Comwell 1989; Birkinshaw & Worrall 1993). The required self-confinement of this 
pressure is easily understood in terms of the magnetically-dominated jet model. The toroidal field 
component B^ provides coUimation of the MHD jet and confinement of high pressure gas via the 
"hoop-stress" {—B^/4nr), which is a part of the magnetic tension force [{B ■ V)B]. 



2.2. Basic Structure of a Current-carrying Jet 

The global picture of a current-carrying jet with a closed current system linking magneto- 
sphere and hot spots, was introduced by Benford (1978) and applied to AGN double radio sources. 
A closed current system includes a pair of current circuits, each containing both a forward elec- 



trie current path (the jet flow itself, with its toroidal magnetic field, toward the lobe), and a return 
electric current path (along some path back to the AGN core). The stability of current-carrying jets 
for two types of return current distributions have been discussed: (1) the returning currents are as- 
sumed to flow within the jet itself (the jets are thermally confined by external medium) (Chiuderi, 
Pietrini, & Torricelli-Ciamponi 1989), and, (2) the returning currents are assumed to flow around 
the jet, such as in a magnetized cocoon (Benford 1978). 



2.3. Stability of HD Jets 

Most of the early work on jet stability concentrated on the purely hydrodynamical (HD) KH 
stability for simple configurations such as "top-hat" velocity profiles (Gill 1965; Hardee 1979, 
1983; Payne & Cohn 1985; Hardee & Norman 1988; Zhao et al. 1992a; Hardee & Stone 1997) 
and other flows (see e.g., Birkinshaw 1991; Ferrari 1998, and references therein). Two types 
of KH waves are considered disruptive in jets: the fundamental surface wave and the reflected 
body wave (Bodo & Ferrari 1982; Zhao et al. 1992a). Surface mode waves are excited in the 
presence of non-zero velocity gradients and/or discontinuities at the boundaries between the body 
of jets and external medium. Body waves propagate through the body of a fluid or only exist in 
the interior of the medium, such as acoustic/magnetosonic waves. The existence of body waves 
does not depend on the presence of boundaries, but the reflection and/or refraction of body waves 
by such boundaries leads to the growth of modes that could eventually cause the disruption of jets 
(Zhaoetal. 1992a; Hardee, Clarke, & Rosen 1997). 

Beginning with the pioneering work by M. L. Norman and his co-investigators (Norman et al. 
1982), numerical simulations have been performed to investigate the nonlinear development of KH 
instabilities for propagating jets (Kossl & Miiller 1988; Cioffi & Blondin 1992; Massaglia, Bodo, 
& Ferrari 1996; Carvalho & O'Dea 2002a,b; Krause 2003). Norman etal. (1982) performed the 
axisymmetric 2-D HD simulations and found that the jet is decelerated by a Mach-disc shock wave 
front that is, in general, much stronger than the bow shock of the jet. Backflow from the working 
surface builds an extensive cocoon or lobe surrounding the jet. Hypersonic jet flow itself is largely 
stable to KH instabilities, which grow slowly compared to the jet propagation speed. The transonic 
lobes and cocoons, however, are generally KH instable and mix with the external gas at a ragged 
boundary. 

Another type of numerical approach investigates the nonlinear development of KH instabili- 
ties in pre-formed ''equilibrium" jets which begin in radial force equilibrium. A thermally confined 
top-hat velocity profiled jet, surrounded by a uniform unmagnetized external medium, was inves- 
tigated with both 2-D (Norman & Hardee 1988; Zhao et al. 1992b; Bodo et al. 1994, 1995; 
Stone, Xu, & Hardee 1997) and 3-D simulations (Hardee & Clarke 1992; Bodo et al. 1998; Xu, 



Hardee, & Stone 2000; Micono et al. 2000). These jets are perturbed at the surface layer between 
jet flow and external medium, and axisymmetric and non-axisymmetric surface and body modes 
of KH instabilities sometimes grow under a variety of conditions. 



2.4. Stability of MHD Jets 

2.4.1. Kelvin-Helmholtz instabilities 

The MHD KH instability of jets has been examined using linear stability analysis (Ray 1981; 
Ferrari, Trussoni, & Zaninetti 1981; Cohn 1983; Fiedler & Jones 1984; Bodo et al. 1989; Hardee 
et al. 1992; Bodo et al. 1996). In general, surface non-axisymmetric modes (m > 0) are stable 
against MHD KH instability during sub-Alfvenic flow. However, in super- Alfvenic but trans-fast 
magnetosonic flow, they can be unstable. The surface symmetric mode (m = 0) is predicted to 
be MHD KH unstable in sub-Alfvenic and super- Alfvenic flow, although with a relatively small 
growth rate (Bodo et al. 1989; Hardee et al. 1992; Hardee, Clarke, & Rosen 1997). However, 
body mode waves can become important and affect the jet interior in the following situations: (1) 
if jet velocity exceeds the super-fast magnetosonic speed (Vfm < ^) or (2) if the flow velocity is 
slightly below the slow magnetosonic speed [CsVa/(C^ + V^)^*^^ < ^j < ^sm] (Hardee & Rosen 
1999). 

The growth of KH instabilities of super-fast magnetosonic jets containing force-free helical 
magnetic field has been shown to be reduced by the presence of a toroidal magnetic field (Appl 
& Camenzind 1992). Similar investigations extended to the sub-fast magnetosonic regime (Appl 
1996) showed that the toroidal field exhibits a destabilizing behavior at small velocities. At least 
for force-free helical magnetic configurations in linear regime, it appears that the KH mode exhibits 
faster growth than CD mode (Appl 1996). Bodo et al. (1996) investigated the axially magnetized 
rotating super-fast magnetosonic jets, and they found that these jets could be stabilized against the 
non-axisymmetric surface mode by jet rotation. If there exists an external wind between the jet 
and the ambient medium, these results are modified considerably. 

Clarke, Norman, & Bums (1986), Lind et al. (1989), and Krause & Camenzind (2001) per- 
formed axisymmetric 2-D simulations of propagating MHD jets with a purely toroidal magnetic 
field and an unmagnetized external medium. Kossl, Miiller, & Hillebrandt (1990) investigated the 
MHD jets with a helical (poloidal -i- toroidal) magnetic field (with the external medium poloidally 
magnetized) using axisymmetric 2-D simulations. The results of these investigations showed that 
a strong toroidal magnetic field makes the backflow weak in the cocoon and, consequently, re- 
duces the growth rate of the KH instabilities. Todo et al. (1992) also performed axisymmetric 
2-D simulations of MHD jets assuming a force-free helical magnetized external medium with 



the injection of super-Alfvenic flows into the computational domain. In their investigation, the 
KH instability is entirely suppressed if the the magnetic field is parallel to the velocity shear be- 
tween cocoon and shroud and if the following criterion is satisfied (e.g., Chandrasekhar 1961): 
Pip2(^i — Vi)^ < (Pi — p2)5^/87r. Here p, V, and B are the density, velocity, and magnetic field 
strength parallel to the velocity, respectively, with subscripts denoting the two different fluids. For 
the present case, subscript 1 corresponds to the cocoon and 2 represents the shroud. Todo et al. 
(1992) derived a qualitative tendency in which the dense external medium (pi/p2 ^ 0. 1) promotes 
the KH instability, whereas the large scale magnetic field suppresses it. 

Recently, using 3-D simulations, Ouyed, Clarke, & Pudritz (2003) investigated the ejection 
and propagation of a MHD jet from a pseudo-Keplerian disk (i.e., the disk is treated as a fixed 
boundary condition which has a Keplerian-like azimuthal velocity profile, but does not allow ac- 
cretion) (see e.g.. Bell & Lucek 1995; Ustyugova et al. 1995; Meier et al. 1997; Ouyed & Pudritz 
1997a,b, 1999). Their results showed that the MHD jet beyond the Alfven surface temporarily 
becomes unstable to the non-axisymmetric (m > 0) KH instability. However, the jet maintains 
long-term stability via a self-limiting process that keeps the average Alfven Mach number within 
the jet to of order unity. This occurs because the poloidal magnetic field becomes concentrated 
along the central axis of the jet, forming a "backbone" in which the Alfven speed is high enough 
to reduce the average jet Alfvenic Mach number to unity. 

Of particular interest are the axially magnetized 2-D slab equilibrium jets (Hardee et al. 1992; 
Hardee & Clarke 1995) and the helically magnetized 3-D equilibrium cylindrical jets investigated 
by P. E. Hardee and his collaborators (Hardee, Clarke, & Rosen 1997). In the super-Alfvenic 
regime, if the jet is also super-fast magnetosonic, then it becomes more stable with increasing 
fast magnetosonic Mach number (Mfm), and the destabilization length (L) varies approximately 
proportional to MpM (L <x MpM^j^ where Rj is the jet radius). Hardee & Rosen (1999) investigated 
helically magnetized 3-D trans-Alfvenic "light" jets (density ratio Pj/pe ~ 0.03, where subscript 
"j" corresponds to the jet itself and "e" corresponds to the external medium). These experienced 
considerable slowing as denser material was entrained following destabilization; provided the jet 
is super-Alfvenic, their KH growth rates also increase as Mfm decreases. However, the jets are 
nearly completely stabilized to these instabilities when the jet is sub-Alfvenic. 

These numerical results show that the super-Alfvenic but trans-fast magnetosonic flow region 
is a potential zone of enhanced KH instabilities just downstream of the Alfven point. They also 
show that magnetic tension can significantly modify the development of the KH instability in non- 
linear regime. An increase of density ratio would stabilized the MHD jets beyond the Alfven point 
and, in general, denser jets have been found to be more robust than their less dense counterparts 
(Rosen et al. 1999; Rosen & Hardee 2000). Hardee & Rosen (2002) confirmed the stabilizing 
influence of a surrounding magnetized wind against the non-axisymmetric KH surface modes of 



the helically magnetized 3-D trans-Alfvenic jets. They concluded that the jets could be stabilized 
entirely to the non-axisymmetric KH surface modes if the velocity shear, AV = Vj — Vg is less than 
a "surface" Alfven speed, Vas = [(pj + pe)(5^5e)/(47rpjpe)]^/^- 



2.4.2. Current- driven instabilities 

Prior to this point in time, much less consideration has been given to CD instabilities than to 
KH instabilities in the investigation of the disruption of the interior of astrophysical jets. This has 
been because, until recently, jets were believed to be in super- Alfvenic or super-fast magnetosonic 
flow {i.e., kinetic energy was expected to exceed magnetic energy). Some analytic studies on 
CD instabilities as a possible explanation for jet disruption had been done (Eichler 1993; Spruit, 
Foglizzo, & Stehle 1997; Begelman 1998; Lyubarskii 1999). Appl, Lery, & Baty (2000) also 
analyzed, for a large range of magnetic pitch {rBz/B(^), the CD linear growth rate for configurations 
with a force-free helical magnetic field and a constant distribution of density and velocity. (These 
assumptions exclude PD and KH instabilities.) They concluded that the properties of the fastest 
growing CD kink (m = 1) mode are nearly independent of the details of the pitch profile. However, 
they also concluded that this was an internal mode that does not cause a significant distortion of 
jet. 

Lery, Baty, & Appl (2000) studied the nonlinear development of CD instability for cold 
super-fast magnetosonic equilibrium jets based on their linear analysis (Appl, Lery, & Baty 2000). 
It was found that the current density is redistributed within the inner part of the jet radius on 
a characteristic time scale of order the Alfven crossing time in the jet frame. Nothing in their 
numerical results indicated a possible disruption of jet by the CD sausage (m = 0) or kink (m = 1) 
mode. However, for numerical reasons, their simulations were limited to early nonlinear phases; a 
full investigation of the nonlinear development of the instability was not carried out. 

Todo et al. (1993) performed 3-D MHD simulations of a model for HH objects in which 
super- Alfvenic propagating jets are injected into a pre-formed force-free helically magnetized am- 
bient medium. They found that YSO jets can be disrupted into a large scale wiggled structure by 
the CD kink instability. In a similar study, Nakamura, Uchida, & Hirose (2001) performed 3-D 
simulations of propagating MHD jets to investigate the formation of wiggled structures in AGN 
jets. They found that the propagation of nonlinear torsional Alfven waves (TAWs) can produce a 
slender jet shape by the "sweeping magnetic twist" mechanism of Uchida & Shibata (Shibata & 
Uchida 1985, 1986; Uchida & Shibata 1985, 1986). In addition, wiggles in the jet can be pro- 
duced by the CD kink instability though the interaction between the TAW and the ambient medium. 
Recently, Baty & Keppens (2002) confirmed the interaction of KH and CD instabilities in heli- 
cally magnetized super-fast magnetosonic equilibrium jets by performing 3-D MHD simulations. 



This nonlinear interaction can contribute to jet survival, and the large-scale magnetic deformations 
associated with CD mode development can effectively saturate KH surface vortices and prevent jet 
disruption. 



2.5. Theoretical Arguments for Poynting Flux Dominated Jets 

There has been a growing recognition in recent years that Poynting flux dominated (PFD) flow 
plays an important role in the jets in AGNs, QSOs, winds from pulsars, and possibly y-ray burst 
sources. This is in strong contrast to the kinetic energy flux dominated (KFD) jets that had been 
expected in these sources, and still are believed to be relevant for YSO jets. The energy carrier of 
PFD jets is primarily the electromagnetic field, while in KFD jets it is the kinetic flux. 

The concept of PFD jets is based on the theory of magnetically driven outflow, proposed (in 
the electromagnetic regime) by Blandford (1976) and Lovelace (1976) and subsequently applied 
to rotating black holes (Blandford & Znajek 1977) and to magnetized accretion disks (Bland- 
ford & Payne 1982). By definition, these outflows initially are dominated by electromagnetic 
forces close to the central engine. In these and subsequent models of magnetically driven out- 
flows (jets/winds), the plasma velocity passes successively through the hydrodynamic (HD) sonic, 
slow-magnetosonic, Alfvenic, and fast-magnetosonic critical points (surfaces). In very strongly 
magnetized flows (Cg ^ Va). the hydrodynamic and slow-magnetosonic points almost coincide, 
and the Alfven and fast points also occur close to each other. The distance between these two clus- 
ters of critical points (HD/slow and Alfven/fast), and therefore the nature of much of the outflow, 
depends on the relative dominance of the advected magnetic field, i.e., on how long the high level 
of Poynting flux domination can be maintained as the flow propagates. 

In many early models of steady, axisymmetric MHD outflow, it was expected that the flow 
would smoothly pass the Alfvenic point (surface) at a distance very close to the central object. 
After that the flow would become super-Alfvenic and, therefore, kinetic energy flux dominated. 
However, according to recent theoretical studies, this picture is probably not correct. Of particular 
interest here is the current-carrying core of the outflow near the central axis (the jet). In steady 
non-relativistic solutions (Fendt & Camenzind 1996; Krasnopolsky, Li, & Blandford 1999) the 
innermost part of the outflow remains sub- to trans-Alfvenic (Ma < 1), i.e., Poynting flux domi- 
nated. Furthermore, in relativistic steady models of MHD outflows from disks, the flow can remain 
PFD for many hundreds to thousands of outer disk radii — tens of parsecs in AGN (Vlahakis & 
Konigl 2001, 2003a,b, 2004). Even in dynamical simulations, a self-organized process occurs 
whereby the Alfven speed increases due to the concentration of poloidal magnetic flux Bp along 
the central axis of the outflow (Ouyed, Clarke, & Pudritz 2003). This reduces Ma near the axis 
and forestalls the passage of the flow through the Alfven point. 
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So, it appears likely that the axial part of the outflow (the current-carrying jet) remains in sub- 
to trans-Alfvenically PFD flow for a long distance. In the case of AGN jets this translates into a 
few to several tens of parsecs. Jets observed by VLBI, therefore, may be Poynting flux dominated, 
making our simulations herein directly applicable. 

Further theoretical development of PFD jets from magnetized accretion disks has been per- 
formed by Begelman, Blandford, & Rees (1984); Lovelace, Wang, & Sulkanen (1987); Li, Begel- 
man, & Chiueh (1992); Lynden-Bell (1996); Romanova & Lovelace (1997); Levinson (1998); 
Colgate & Li (1999); Li et al. (2001); Lovelace et al. (2002); Hey vaerts & Norman (2003a,b,c); 
Lovelace & Romanova (2003); and Vlahakis (2003). Two-dimensional axisymmetric ("2.5-D") 
non-relativistic MHD simulations of PFD jets have been performed for opening magnetic loops 
threading a Keplerian disk (Romanova et al. 1998; Ustyugova et al. 2000). And relativistic 2.5-D 
simulations of propagating PFD jets carrying a toroidal field component only were performed by 
Komissarov (1999) in order to compare their simulations with the non-relativistic simulations of 
Lind et al. (1989). Finally, Li (2000), Tomimatsu, Matsuoka, & Takahashi (2001), and Wang 
et al. (2004) have considered analytically the CD instability (the so-called "screw instability") of 
black hole magnetospheres. 

To our knowledge, however, no study has been made of the nonlinear behavior of PFD jets 
and the related CD instabilities in full 3-D dynamics. 



3. NUMERICAL METHODS 

3.1. Basic Astrophysical Model 

In this paper we study the structure, dynamics, and stability of propagating MHD jets. Our 
simulations concentrate on the region in an AGN in which a coUimated jet-like flow (at much 
greater than the escape velocity) has been established, but that flow still is dominated by magnetic 
forces and has not yet achieved a super-(fast magneto) sonic velocity. The energy carried by the 
jet, therefore, is dominated by Poynting flux rather than kinetic energy flux, and the Alfven speed 
in the flow is much higher than the local sound speed, Cs ^ Va- 

Numerous theoretical investigations of the central engine itself have shown that a rotating, 
magnetic structure can be created a few Schwarzschild radii (10^^^^^ cm) from the black hole 
(Blandford & Znajek 1977; Blandford & Payne 1982; Koide et al. 2002). While the physical 
connection between this central region and the sub-parsec region is still poorly understood, it is 
reasonable to suppose that both strong magnetic and rotational fields will be influential at the base 
of the sub-parsec flow. Our basic jet model, therefore, contains a plasma with a strong, rotating 
poloidal magnetic field at its base. Nonlinear TAWs propagate out along this field, creating and 
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propagating a magnetically driven coUimated outflow via the sweeping magnetic twist mechanism 
of Uchida & Shibata (1986). 

In contrast with the many numerical stabilities studies, in which the magnetized jets are as- 
sumed to be confined thermally by a non-magnetized external ambient medium, we assume a large 
scale poloidal magnetic field in the ambient medium surrounding the jet. The origin of such a 
galactic magnetic field is not yet is fully understood, but its existence is suggested by both syn- 
chrotron emission and Faraday rotation observations. The magnetic field assumed here might 
either be a part of the primordial inter-stellar field (Kulsrud & Andersen 1992) brought into the 
central part of the protogalaxy in the process, or the central part of the amplified field by a galactic 
turbulent dynamo process, which are argued by many authors (for reviews, see Kronberg 1994; 
Han & Wielebinski 2002, and references therein), or a field structure carried out from the central 
engine by a lower- velocity magnetized disk wind. 

Models of the central parsec and sub-parsec regions of galaxies indicate that the ambient 
medium is characterized by strongly decreasing gradients in density and pressure and, presumably, 
magnetic field. We therefore are especially interested in the behavior of PFD jets in decreasing 
density atmospheres. 



3.2. MHD Equations and Simulation Code 



In the present study we assume non-relativistic ideal MHD and neglect the effect of gravity. 
We solve the nonlinear system of time-dependent MHD equations numerically in 3-D Cartesian 
coordinate system (x, j, z): 



| + V.(pV)=0, 



f + (v.v)v' 



-Vp + JxB, 



dp 






(1) 

(2) 

(3) 
(4) 



Here, p is the mass density, p the gas pressure, V the fluid velocity. B is the magnetic field and 
7 = (V X B)/4n is the corresponding current density. F is the ratio of specific heats (a value of 5/3 
is used). 

We normalize all physical quantities with the unit length scale Lq, typical density po, typical 
velocity Vq in the system, and their combinations, e.g., p' = p/po, etc. The normalizing factors are 
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shown in Table 1. A factor of 4n has been absorbed into the scaling for both the magnetic field 
B and the current density J. This normalization does not change the form of the basic equations 
(1) - (4), and hereafter, we will use the normalized variables and will omit the primes on physical 
quantities. (We will write the dimensional variables with subscript "0": po, for example.) The 
system of dimensionless equations is integrated in time by using a two-step Lax-Wendroff scheme 
(Rubin & Burstein 1967) with artificial viscosity term (Lapidus 1967). The original MHD code 
was developed by K. Shibata and extended by his co-workers (e.g., Shibata 1983; Shibata & 
Uchida 1985, 1986; Matsumoto et al. 1996). Parallelization of the code was done using MPI. 

The total computational domain is taken to be \x\ < Xm^x, \y\ < Jmax, and Zmin < z < Zmax, 
where Xmax, Jmax — 16.0, Zmin — —1-0, and Zmax — 20.0. The numbers of grid points in the sim- 
ulations reported here are A^x x A/j, x A^^ = 261 x 261 x 729, where the grid points are distributed 
non-uniformly in the x, y, and z directions. The grid spacing in the initial propagation region near 
the jet axis is uniform, Ax = Ay = Az = 0.025 for |;c|, \y\ < 1.5 and |z| < 16.0 (121 x 121 x 681 
cells are assigned here), and then stretched by 5% per each grid step for the regions \x\, \y\ > 1.5 
and |z| > 16.0. Forty uniform cells are distributed across the initial jet diameter in the transverse 
direction, L ~ 1 at z = 0. Low-resolution exploratory computations were performed on the Jet 
Propulsion Laboratory Origin 2000 machines. High-resolution 3-D computations were performed 
on the FUJITSU VPP 5000/32R at the National Astronomical Observatory in Japan (9.6 GFLOP/s 
peak speed per node) requiring about 6 CPU hours each. 



3.3. Initial Conditions: Non-uniform magnetized atmospheres 

In the non-uniform, magnetized atmosphere we adopt a current- (and therefore force-) free 
magnetic configuration (7 = 0), by placing one pair of current loops on both the upper and lower 
z-boundaries of the computational domain. This constrains the magnetic field to have only a field z- 
component at z = Zmax, thereby avoiding numerical errors when physical MHD waves pass though 
the upper (z = Zmax) boundary. Such a field configuration does not disturb the force equilibrium 
of the hydrostatically stable atmosphere because J x B = 0. The field configuration chosen is 
explicitly given in the Appendix. 

For our initial density distribution, we assume that p varies as a power of the magnetic field 
strength 

poc|B|«, (5) 

where ais a free parameter in this paper. If a = 2, the Alfven speed Va(= I-^I/a/P) will be constant 
throughout the computational domain. If a 7^ 2, the Alfven speed will decrease (a < 2) or increase 
(a > 2) with distance from origin (0, 0, 0). This power law model embraces several accretion and 
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collapse models for the formation of the interstellar medium (ISM). 

• For example, if the ISM had formed by the conservative spherical contraction of a pro- 
togalactic gas cloud, the magnetic field strength |fi| would be amplified as |fi| oc p^/^ or 
a = 3/2. Combining this value for a with our initial magnetic field distribution, the po- 
lar (axial) distribution of Alfven velocity Vaz{z){= B^/ ^) would gradually decrease as 

VUz) - z-^'^ [B,{z) - z-^, p(z) oc z-\ 

• In the case of interstellar magnetized clouds (isothermal gravitational contractions), 2 < 
a < 3 is theoretically (Mouschovias 1976), numerically (Scott & Black 1980), and obser- 
vationally inferred (Crutcher 1999) rather than a = 3/2. However, a varies widely with 
position (r,z) in the contracting gas, and is smaller than 2 near the central (z) axis, because 
the field there is being "dragged along" by the collapse in the r-direction. This leads to a 
large amplification of the field with no corresponding large increase in p; thus a becomes 
small (Mouschovias 1976; Scott & Black 1980). a ~ 2 [p(r) oc 5^(r)^] on the equatorial 
plane and a ~ 1 [p(z) oc B^{z)] on the central axis of collapsing gas, also have been found 
with numerical simulations (Tomisaka 1996). 

• A third possibility is that the ISM may have formed in an advection-dominated accretion flow 
(Narayan, Mahadevan, & Quataert 1998), where p °^ r^'^l^ and \B\ ^ r^^/^, or a = 6/5. 

We therefore choose the following two representative cases: ''decreasing Va" (cx = 1) and 
''constant Va" (oc = 2) as the initial ambient medium. 

For the initial gas pressure distribution we make the polytropic assumption 

P °- P^, (6) 

where F is the polytropic index (we use F = 5/3 throughout this paper). The sound speed Cs 
(= a/Fp/p oc r^/^, with T being the temperature) will decrease with distance from origin. The 
atmosphere is artificially confined to prevent it from expanding under its own pressure gradient 
by imposing a pseudo-gravitational potential (^) designed to 'hold on to' the atmosphere without 
significantly impeding the advancing jet (Clarke, Harris, & Carilli 1997): 

Here, p and p are the initial gas pressure and density distributions. By design, the quantity p^ 
exactly cancels the initial gas pressure gradients in the stratified atmosphere, V(p — p^) = 0. As- 
suming a current-free (J = 0) field, the right-hand side of equation (2) exactly equals zero. 
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We now consider the "plasma-(3" parameter: 

M^=4 (8) 

^ |fl|2 ry2 ^ 

as the ratio of the gas to the magnetic pressure. In all the situations we investigate in this paper, the 
thermal energy density is much less than the magnetic energy density, so |3 ~ C^/V^ ^ 1 in the 
whole computational domain. 



3.4. Boundary Conditions 

In the lower "boundary zone" [zmin < z < 0, where r = (jc^ +y^)^/^] we set the velocity field 
to be 

V^=V^{r,z)^ + V,{r,z)z. (9) 

This zone is the injection region, in which the hypersonic but sub-Alfvenic flow (i.e., the PFD jet) 
is formed. Equation (9) represents a cylindrical MHD jet, powered by a nonlinear TAW, entering 
the upper region (z > 0) of the computational domain. In all our simulations the injection speed of 
the PFD jet in the boundary zone is sub-Alfvenic, but super-slow magnetosonic. 

To stabilize the numerical calculations in the injection region, we employed a "damping zone" 
to prevent unnecessary reflections and interactions from below. The disturbances in all physical 
quantities except for the magnetic field are damped at each time steps as follows: 

Q{x, J, z, ?) = { 1 - fQ{z)}Qix, y, z, t) 

+fQ{z)Q{x,y,z.O), (10) 

The damping factor /q for quantity Q is defined such that it is equal to zero at the upper end of the 
damping zone, z > Z2, and to unity at the lower end, z< Zi- We set the values to zi = —0.975 and 
Z2 = —0.025, respectively. We use about 40 axial cells in this boundary zone. Note that we leave 
the magnetic field B{xj y, z, t) unmodified so that the induction equation (3) still holds. 

At the bottom of this zone (z = Zmin), we impose symmetry for p, Vx, Vy, B^, and p, and 
antisymmetry for V^, B^, By. The free time-dependent evolution of the flow physical variables 
occurs only in the upper part of the computational domain (z > 0). 

At the outer boundaries (x = Xmin, x = Xmax, y = ymin, y = Jmax, and z = Zmax) we set free 
conditions for each quantity Q as follows Shibata (1983): 

dAQ _ dAQ _^^Q _Q 
dx dy dz 
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^Q = Q{x,y,z,t+M)-Q{x,y,z,t). (12) 

4. RESULTS AND DISCUSSION 

In this section we provide a detailed description of the evolution of four different simulations 
of PFD jets. We will concentrate on the solutions in the "jet-propagation" region, —2.25 <x,y< 
2.25, and 0.0 < z < 1 8.0. An important unit time scale in these dimensionless systems is the Alfven 
crossing time Xa (= L/Va), normalized with Xao = Lq/Vaq- 

4.1. The Four Initial Models 

We have carried out a number of simulations with varying physical conditions, evolving four 
distinct models: shallow-atmosphere models (A) and steep-atmosphere models (B), each with a 
highly and a mildly Poynting flux dominated flow case. The values of their physical parameters 
listed in Table 2 and discussed more fully below. 

The main purpose of this investigation is to study the effects of a decreasing density, mag- 
netized ambient medium on the growth of current-driven instabilities in dynamically propagating 
PFD jets. Figure 1 shows the initial {t = 0.0) distribution of the physical quantities along the z-axis. 
For the A models (a = 1; decreasing Va), p and B^ decrease gradually along the z-axis and asymp- 
tote to ~ z ^ for z > 1. On the other hand, for the B models, B^ decreases in the same manner as 
model A, but p decreases faster than ~ z^^ for large z (cx = 2; constant Va)- We set the plasma 
beta (P = 10^^) at the origin for all of the models. The ratio Va/Q decreases slightly along z-axis 
for model A but gradually increases for model B. 

We also consider two different types of flow in each atmosphere model, i.e. the ratio of 
Poynting flux F^xs to the total energy flux Ftot in MHD jets: 

1. Highly PFD jets (Fe^b/Fioi ~ 0.9) (Case 1) 

2. Mildly PFD jets (F^xs/Ftot ~ 0.6) (Case 2) 

Figure 2 shows the time variation of the ratio of the injected energy fluxes into the upper "evolved 
region" from the lower boundary zone. Note that the quasi- stationary PFD flow is injected through- 
out the time evolution of the system. The net amount of injected Poynting flux for the mildly PFD 
jets is almost equal to the case of the highly PFD jets. However, the amount of kinetic energy flux 
injected in Case 2 is considerably larger than in Case 1. 
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4.2. Comparative Overview of Simulation Evolution 

Before considering our numerical results in detail, it is instructive to give a brief comparative 
overview of the time development of the jet systems in the four models. This can be done by 
inspecting the logarithmic density p and the velocity field (Vx, V^) normalized by Vao- Figures 3 - 
6 display these quantities at various times in the simulation evolution using two-dimensional x — z 
slices at y = 0. Note that these slices contain the jet (z) axis. 



4.2.1. Early jet evolution: propagation of the nonlinear MHD waves 

At early times in the evolution {Top part of each Figs. 3-6), the well-coUimated helically 
PFD jets powered by the TAWs advance into the poloidally magnetized ambient medium. Figure 
7 displays for all models the axial distributions of the Alfven and sound speeds, as well as each 
component of the jet velocity, in the x — z plane close to the jet axis. The propagation of three 
distinct fronts of compressive MHD waves is clearly visible for all models: a fast- mode MHD 
wave (at the head of the PFD jet or TAW) and a pair of slow-mode MHD waves further upstream. 
We characterize these waves as propagating in both the forward (F) and reverse (R) directions 
relative to a reference frame that co-moves with the jet, and give them the names forward fast- 
mode (F— F), forward slow-mode (F— S), and reverse slow-mode (R— S) waves. In addition, a 
contact discontinuity travels in this co-moving frame between the two slow-mode waves. Because 
the jet is injected at a super-slow magnetosonic velocity in the boundary zone, the forward and 
backward (reverse) propagating slow-mode compressive waves seen in each of the models quickly 
become slow-mode shock waves early in the simulation. 

In models A - 1 and A - 2, due to a gradual decreasing ambient Va, the F— F compressive wave 
front decelerates and its wave amplitude gradually increases as it propagates forward. Through this 
nonlinear process, this front steepens into di fast-mode MHD bow shock when the phase velocity 
becomes super-fast magnetosonic (see also the 2nd panels in Fig. 3 and 4). 

On the other hand, only a very low amplitude F— F compressive wave front can be (barely) 
seen in models B - 1 and B - 2 in Figure 7. This wave front never reaches a super- fast magnetosonic 
velocity during our simulations of models B. The velocity amplitude associated with the F— F 
compressive wave becomes much weaker as the ratio Va/Q in the ambient medium [= (2/r(3) ^' ^] 
along z-axis increases. The phase velocity of a fast-mode compressive MHD wave is Vph = VpM ~ 
Va when Cl <^ V^. Because the initial distribution of the Alfven velocity in Models B - 1 and B - 
2 is constant with z (Va = 1-0), the F— F compressive wave will propagate quickly compared with 
models A - 1 and A - 2. In the B models the wave reaches z = 18.0 at a time t ~ 18.0. 

Several features of PFD MHD jet flow have counterparts in the well-known supersonic hydro- 
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dynamic jet flow. The F— F compressive/shock wave, of course, propagates through the external 
medium at a speed faster than the other two shock waves, playing the same role as the bow shock 
does in HD jets. The F— S shock wave propagates at a slower speed, further compressing and heat- 
ing the ambient medium. This second shock wave has no direct counterpart in HD jets, although it 
acts similarly to the bow shock, especially when the F— F front is a compressive wave. The reverse 
(R— S) shock wave plays the same role as the Mach disk does in HD jets, decelerating and heating 
the magnetized jet flow itself. From the point of view of the jet material, the flow can reach the 
R— S shock wave and enter the compressed region, but it can never reach the F— S shock wave. 
Similarly, the external ambient material crossing the F— S shock wave also reaches this heated re- 
gion, but can never reach the R— S shock wave. Instead, ambient material is entrained between the 
two slow-mode shock waves, along with the accumulated jet material. The contact discontinuity 
between the two accumulations, therefore, plays the same role as the contact discontinuity in HD 
jets. We often identify this interface as defining the rest frame of the jet flow. 



4.2.2. Intermediate jet evolution: energy conversion though the MHD shocks 

In the "intermediate stages" of the flow {2nd panel each in Figs. 3-6), the gas is weakly 
compressed in the z direction across the first (F— F) compressive/shock wave, and then is strongly 
compressed further across the second (F— S) shock wave. In addition, the gas crossing the third (R— 
S) shock wave also is expanded by a large ratio, but in the opposite direction. Because the greatest 
compression occurs when the gas crosses the slow shock waves, the material undergoes extensive 
heating and accumulates in the region between these two wave fronts. In the shallow-atmosphere 
models A the structure of the pair of slow shock waves does not change during the simulation. 
However, in models B this structure undergoes an expansion in the transverse direction because of 
the steeply decreasing gradient in the external ambient pressure. 

Figures 8 and 9 show snapshots of Vy and By (similar to Fig. 3 - 6) at the time when the 
head of the PFD jet reaches z = 10.0 in each of the models. (Vy and By correspond to V^s^ and B^ 
in cylindrical coordinates.) Below the flow through each compressive/shock wave is described in 
the shock frame (in the negative z direction for the F— F and F— S shock waves and in the positive 
z direction for the R— S shock wave). We now can compare the energy conversion in each of Figs. 
8 and 9: 

1. The toroidal component of the magnetic field {By) increases strongly as the flow crosses 
the F— F shock wave front (in the -z direction; models A - 1 and A - 2). However, that 
component increases only weakly across the F— F compressive wave (models B - 1 and B - 
2). Therefore, the toroidal (rotational) kinetic energy of a PFD jet is strongly converted into 
toroidal magnetic energy across the F— F shock wave in the shallow-atmosphere models A, 
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but not so much in the steep-atmosphere models B. 

2. The toroidal component of the magnetic field decreases across the F— S shock wave, and 
the toroidal component of the velocity field also decreases across it. Therefore, the toroidal 
kinetic and magnetic energy of a PFDjet are converted into thermal and axial kinetic energy 
across the F— S shock wave. 

3. The toroidal component of the magnetic field also decreases as the flow crosses the R— S 
shock wave (in the +z direction). However, unlike the first slow-mode shock wave, the 
toroidal component of the velocity field does not decrease across the second. It actually 
increases as material flows into the region between the two slow-mode shock waves. (Note 
the relative strength of the toroidal component of the velocity beyond the F— S shock wave 
and behind the R— S shock wave in the frame co-moving with the contact discontinuity for 
Case A - 2 in Fig. 8.) Therefore, axial kinetic and toroidal magnetic energy are converted 
into thermal and toroidal kinetic energy across the R— S shock wave. ^ 

So, we expect that the magnetic field between the second and the third shock waves will be 
less twisted, whereas that between the first and second, and behind the third, shock waves will be 
more twisted. This physical picture is similar to the results of 1.5-dimensional MHD simulations 
performed by Uchida et al. (1992). 



4.2.3. Nature ofPFDjets as current-carrying jets 

We now consider the behavior of PFD jets as current-carrying systems. Figure 10 shows 
snapshots of the axial current density J^ (similar to Fig. 8 and 9) at the time when the head of the 
PFD jet reaches z = 10.0 in each models. Clearly, the flow displays a closed circulating current 
system, in which one path occurs close to the central axis, and co-moves with the PFD jet (the 
'forward jet current density" J^^) and another conically-shaped path that flows outside (the ''return 
current density" J'^^). 



'At first it would appear that the behavior of two slow-mode shock waves is strangely asymmetric. However, if 
one transforms to the translating and rotating frame of the contact discontinuity between these two shock waves, one 
finds that the behavior of all energy components (axial and toroidal kinetic, toroidal magnetic, and thermal) are all 
precisely symmetric. In particular, at both slow-mode shock waves toroidal kinetic energy decreases in this frame, and 
is converted into thermal energy, as the flow crosses into the region between them. There is only one real asymmetry: 
the change in the sense of rotation across the two slow-mode shock waves is anfi-symmetric. This antisymmetry is 
cause by the axial speeds of the two shock waves being opposite in this frame while the handedness of the magnetic 
field pitch, and therefore V,^, is the same. 
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Figure 11 shows (at the same times as in Fig. 10) the transverse profile of the force-free 
parameter Xff (= J ■ B/\J\ \B\) and the axial current density J^ near the lower (z = 2.2) and upper 
(z = 8.3) regions. The diagram shows that the jets relax into radial force balance with the ambient 
medium: there is a current-carrying core J^^ in the jet itself and a return-current skin J'^'^ outside, 
both of which are in approximate force-free equilibrium (7 x B ~ 0). (The reversal of sign from 
?iff = — 1 to 1 indicates a change in sign in the direction of axial current flow.) Moreover, an almost 
(axially) current-free (7^ — 0) sheath lies between them, as can be seen in panel (2-2) of Fig. 11 
(0.4 < \x\ < 1.0). The net axial current 4 must be zero to avoid the accumulation of electrical 
charge at the end point: 

f ''max 

k= 2'!l{Jf+Jl^)rdr = 0. (13) 

/"max is located at the radius where B^{rmw) = 0. This is similar to the transverse analytical struc- 
ture used by Lind et al. (1989) as inflow conditions for their axisymmetric, toroidal field-only 
simulations. At intermediate times, therefore, our three-dimensional numerical simulations have 
the same dynamical structure as a current-carrying equilibrium jet in force balance. 

A remarkable feature of the Jf distribution is the high axial current density in the region 
5 < z < 10 in models A - 1 and A - 2 {Top and 2nd panels in Fig. 10), and low axial current 
density there in models B - 1 and B - 2 {3rd and Bottom panels in Fig. 10) [see also (2-2) of 
Fig. 11]. This is directly related to how the toroidal component of the field B^ increases across the 
first (F— F) shock/compressive and third (R— S) shock waves {or decreases across the second (F— S) 
shock wave]. A strong accumulation of 5(|) produces a large hoop-stress (— 5?/r), pinching the jet 
towards the central axis. As a result, the gradient of the B^^f in the transverse direction becomes 
larger, and jf along the central axis is therefore enhanced. 

As discussed above, we identify an almost (axially) current-free region between the forward 
jet current J^^ and the return current J'^'^ densities, where J^ is close to [see Fig. 10 and (2-2) of 
Fig. 11]. In this region the transverse component of the field By and velocity Vy have a maximum 
value and the magnetic lines of force are highly twisted (with the ratio By/B^ much larger than 
that in the central jet). We will refer to this region as the "external magnetized wind". (The total 
outflow consists of both the jet and the wind. Here we distinguish the wind from the jet by the 
fact that it is not current-carrying.) In a later section we will discuss the wind's effects on the 
stabilization of the jet against disruption by KH instabilities. 



4.2.4. Final evolutionary stages: destabilization of PFD jets 

In the late stages of jet evolution {3rd and Bottom panels of each of Figs. 3-6), the jets 
are deformed into wiggled structures in both the density and velocity fields. For the shallow- 
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atmosphere, highly PFD model A - 1, the jet is distorted both beyond the second shock wave and 
behind the third shock wave. In addition, as the jet instabilities become fully developed, the high 
temperature region containing the pair of two slow-mode shock waves and the contact discontinuity 
between them, is fully disrupted also. 

On the other hand, for models A - 2 (the shallow-atmosphere, mildly PFD case) and B - 1 
(the steep-atmosphere, highly PFD case) the jets are distorted behind the second shock wave only. 
The high temperature region between the two slow-mode shock waves appears unaffected. 

The instability appears to have the slowest growth rate in the steep-atmosphere, mildly PFD 
case (model B - 2), with only a slight distortion appearing behind the R— S shock wave. 

Figure 12 shows snapshots of J^ in the "final stages" of all the models. Here we also see that 
the jet current density jf is distorted into the wiggled structure seen in the density and velocity 
fields of models A - 1, A - 2, and B - 1. In Fig. 13 we also show (in the final stages in all models 
and in the x — z plane close to the jet axis) the distributions of the Alfven and sound speeds, as well 
as each component of the jet velocity. We identify the following features in each of the models: 

• The axial speed of the jet V^ in model A - 1 is still sub-Alfvenic except at the jet head itself 
(the head is a F— F shock wave, and therefore, super- Alfvenic and super-fast magnetosonic); 
that is, the highly PFD jet remains Poynting flux dominated in this shallow atmosphere. 

• V^ in model A - 2 is super- Alfvenic beyond z > 2 (the jet head is also a F— F shock wave); 
that is, the mildly PFD jet has switched to a mildly kinetic energy flux dominated (KFD) jet. 

• V^ in model B - 1 is trans-Alfvenic beyond z > 10, where the jet is distorted; that is, the 
highly PFD jet in this steep atmosphere has switched to an equipartition state between the 
kinetic and Poynting fluxes. 

• Vz in model B - 2 is super- Alfvenic beyond z > 2; that is, the mildly PFD jet in the steep 
atmosphere has switched to a KFD jet. 

We conclude, therefore, that jets propagating in the trans-Alfvenic region, before they become 
kinetic energy flux dominated (Va ~ ^fm ^Vj)^ can be deformed into wiggled structures. We 
consider the physical mechanism for these destabilization in the following sections. 
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4.3. Nonlinear Growth of CD Instabilities 

4.3.1. Computation of the power spectrum 



To better understand the nonlinear behavior of CD instabilities in our jet simulations, we begin 
by defining the forward jet current density J^^ to be equal to the total current density J in those 
places where J^ has a negative sign and zero otherwise: 



7J^=/ (^^<'^ (14) 

"^.0 (7,>0)- ^'^^ 



We distinguish J^'^ from the jet return current density 

J (Jz>0) 



(7^<0) ' ^^^^ 

which differs in sign from J^'^^. Note that /J'^ (and 7'^'^) is a vector and can have additional radial and 
azimuthal components JI^ and 71*^ with any sign. (These components are computed by transforming 
to a cylindrical coordinate system from the simulation Cartesian coordinates.) 

We now proceed to analyze the modal structure of 7-''^ by forming the volume-averaged Fourier 
transform of its magnitude 

Ji^{m, k) = — [[[ l/J^le'M+fe) rdrdi^dz, (16) 

Vci JJJv^i 

1 /9 

where |7J'^| = [{JI^)^ + (J^^)^ + (Jz^)^] ^nd Vci is the annular volume that encloses 7-''^ 

fb f2n rzh 

Va= / / rdrd(^dz. (17) 

Jr.^ Jo Jza 

We use the values r^ = 0.075, rb = 2.25, Za = 0.0, and Zb = 18.0. 

J^^{m, k) is now a function of the azimuthal mode m and axial wave number k = In/X (corre- 
sponding to a characteristic wave length X). Like 7^'^, 7^^^ is a function of time. 

Finally, we identify the power spectrum as the squared amplitude of 7^^^ 

|7J'=(m,/t)|2=|Re[7J'^(m,/t)]} + {lm[7J'=(m, /t)] | , (18) 

Re[7J'^(m,/t)] 

= —111 \J^''\cos{m(^ + kz)rdrd(^dz, (19) 

Vci JJJv.i 

Im[7J'^(m,/t) 



— /// \ji^\sm{m(^ + kz)rdrd(^dz. (20) 

Vci JJJv^i 



22- 



For each model we now will examine the time-dependent behavior of the power spectrum and 
investigate the nonlinear growth of each mode. 



4.3.2. Growth of azimuthal modes 

We first will consider only azimuthal modes (m) — i.e., the power spectrum \J^^{m)\ for 
infinitesimal wave numbers (A: ^ 0). Several modes (m = 1 — 4) are plotted as functions of time in 
Fig. 14. Note that, while we have not imposed any explicit perturbations to our jets, they are subject 
to the implicit perturbations created by the 3-D Cartesian grid. As a result of our resolving an 
initially, rotating cylindrically symmetric object on a Cartesian coordinate, the m = 4 mode appears 
before the physical growth of the other asymmetric modes. We find, however, that the power in 
the m = 4 mode remains constant or decreases with time on all models. We therefore identify the 
level of power in the m = A mode as our (conservative) noise level and identify all power greater 
than that (and in any other mode) as spectral power generated by the physical evolution of the jet. 

In all models the kink (m = 1) mode grows faster than the other higher order modes (m = 2, 3). 
In Fig. 14 the time at which the kink mode appears above the "noise" level in each of the models 
is as follows: ? ~ 18 (A - 1), ? ~ 23 (A - 2), ? ~ 23 (B - 1), and r ~ 17 (B - 2). In addition, in A - 
1 both the m = 1 and m = 2 modes appear around t ~ 30, and the m = 3 mode saturates at or below 
the noise level. In A - 2, the m = 2 and 3 modes track one another and saturate at our chosen noise 
level. In B - 1 and B - 2, the m = 2 and 3 modes also track one another and saturate at a level 
lower than that of the m = A mode. 

The m = 1 mode for each model exhibits approximately exponential growth on a dynamical 
time scale. A linear fit using a least-squares method has been performed on the time- sequenced 
data to extract a linear growth rate 

^ , , Jln|7>(m)|2 

Im(a))~ ' ^ ' . (21) 

The estimated growth rates are given in Table 3 and shown on Fig. 14. The growth rate appears 
to be strongest for A - 2 [Im ((o) ~ 0.54], and weakest for A - 1 [Im((o) ~ 0.31]. As seen in Fig. 
14, the instability appears to have saturated in the final stages for both A - 2 and B - 2 (the mildly 
PFD models), but still appears to be growing for both of the highly PFD models A - 1 and B - 1 . 

How do these growth rates compare with the Alfven crossing time scale, which is typically 
the time scale on which pure CD modes can develop (Begelman 1998; Appl, Lery, & Baty 2000; 
Lery & Frank 2000)? The spatial length scale of the CD instability in our results (in the radial 
direction) is of the order of a jet width, L ~ 1.0 (see Fig. 12). The local Va in the distorted jet 
portion is of order 0.3 — 0.8, as seen in Figs. 7 and 13. So, the inverse of the Alfven crossing time 
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x^ in the distorted jets is approximately given by 



-1 Vp^ 
L 



Xa =-r^~0.3-0.8. (22) 



This is, in general, consistent with the time scales of the growing m = 1 modes derived from Fig. 

14, 

Im(co)^0(xX^). (23) 



4.3.3. Development of the m=\ helical kink 

Figure 15 shows a space-time (z, t) plot of the ratio of the transverse and axial magnetic field 
components —By/B^ for model A - 1, illustrating the dynamical evolution of the helically kinked 
magnetic field. Propagation of the jet head, corresponding to a F— F compressive wave front, can 
be seen. The wave gradually decelerates due to the decreasing ambient Va until ? ~ 15; after 
that it propagates with an almost constant phase velocity. The amplitude of the F— F compressive 
wave grows significantly with time and could steepen into a F— F shock wave (bow shock) with 
subsequent evolution. After t ~ 20, the kinked part of the field appears behind the R— S first, and 
later, it appears in the upstream region. Finally, the flow downstream (between the F— F and the 
F— S) is also kinked. We can see that the kinks in these three regions grow independently, and the 
peaks of the kinked magnetic field simply co-move with the jet. These characteristics indicate that 
the CD kink mode does not propagate in the jet rest frame (Appl 1996). This behavior is quite 
different from that of the KH instability. In the latter case, the disturbance propagates backwards 
in the jet rest frame, because the phase velocity of the KH modes is smaller than the jet velocity, 
which must be super- Alfvenic to stimulate the KH unstable modes. 

Figure 16 shows snapshots of the radial specific power generated by the Lorentz force on two- 
dimensional x — z slices at J = 0. This is defined as the work done per unit mass in the r-direction 
in a cylindrical coordinate system (r, (|),z): 

(7xB),K/p, (24) 

where 

{JxB)r = -J,B^+J^B^. 

If we consider only radial variations in the Lorentz force, then 

-7,5^ = -^|-(r5c^), 
r or 
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and 

There is a strong correlation between this power distribution and the wiggled structures seen in 
A - 1, A - 2, and B - 1. This is consistent with our early finding that the disruption occurs only 
locally where the MHD jet is CD unstable. That is, distortions do not propagate. These snapshots 
also lead us to conclude, as we did in Fig. 14, that the CD instability is still growing for the highly 
PFD jets A - 1 and B - 1, but is almost saturated for the mildly PFD flows A - 2 and B - 2. This 
means that the magnetic field in the flow should be in "radial" force-free equilibrium (7 x B ~ 0) 
in the saturated models. 

We conclude that CD instabilities grow during the simulation evolution even without the ad- 
dition of small perturbations to the system. This is because we have performed computations with 
high enough resolution to amplify asymmetric modes from the numerical background noise. We 
close this section by noting that the instabilities reported here do not arise if we choose a grid 
spacing that is two or more times larger than the current one. The effects of numerical diffusion 
appear to damp out the growth of the CD instability in this case. 



4.4. Stabilization of the CD Instability with Rapid Rotation 

4.4.1. The Classical Kruskal-Shafranov Criterion 

According to the well-known Kruskal-Shafranov (K-S) criterion (Shafranov 1957; Kruskal 
et al. 1958), which is based on a skin-current model, a magnetized flux tube (along the z axis) is 
stable to the CD kink (m = 1) mode as long as the magnetic twist angle $(r) is below some critical 
value 4>crit, 

4> r =-^< 4>ent, (25) 

rBz 

where, L is the length of the current-carrying magnetic flux system, r is the radius of the flux tube, 
and B^ and B,^ are the axial and azimuthal field components in cylindrical coordinates, respectively. 
Many theoretical and numerical investigations have been performed on the current-driven stability 
problem of solar coronal loops, and the K-S criterion plays an important role in those studies. 

In the original K-S criterion, 4>crit is set equal to 2n. The effects of line-tying, i.e., anchoring 
the axial boundary conditions on the magnetic field in a high-density photosphere, however, raise 
the stability threshold. The modified critical values for a force-free configuration are between 2k 
and IOtt (Hood & Priest 1979; Einaudi & Van Hoven 1983). 



25 



For a jet configuration, wliicli lias L:^ r,we replace L by a spectrum of axial wavelengths X, 
and we define the critical wavelength to be the one that is marginally stable 

^crit = (26) 



The K-S stability criterion (eq. 25) then becomes 

^<Ut, (27) 

so that wavelengths greater than Xcrit are unstable to the CD instability. From a computational 
point of view, X must lie in the range 

A<X<L{t), (28) 

where A is the minimum grid width and L{t) is time-dependent jet length, with maximum value 
^max(0 dependent on the computational domain. In practice, however, shock waves and other 
MHD phenomena provide effective domain boundaries, so that the maximum value of X could be 
only a fraction of the entire jet length. 



4.4.2. Force Balance in Rotating Jet Systems 

Because real astrophysical jet systems may be rotating, the K-S criterion may not be an 
adequate stability condition for the CD kink mode. Rotation will modify the radial force balance 
inside jet. Although rotation is often neglected in theoretical models of jet equilibrium (Begelman 
1998; Appl, Lery, & Baty 2000), it could play an important role in the stability analysis. Let us 
consider the time-independent, force-balance equation (2) in the radial direction. In an inertial 
frame co-moving with the jet we have 

-VP + yxfi + fpaV = 0, (29) 

where Q. is the equilibrium angular velocity and f is the unit vector in the r direction. It is clear 
that the axial flow velocity V^ does not affect the radial force balance, but the presence of an 
azimuthal velocity component can have a significant influence on the equilibrium, modifying the 
competition between the pressure gradient and the Lorentz force. Moreover, our PFD jets are 
strongly magnetized ((3 ^ 1), so the pressure gradient term on the left-hand side of equation (29) 
can be neglected. This yields 

* -^ + ^ = 0, (30) 
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where Vt^ is the equilibrium azimuthal velocity. The first term on the right-hand side in equation 
(30) is the magnetic pressure gradient in both the azimuthal and axial components. These forces 
are positive or negative depending on whether B^ and B(^ decrease or increase with r, respectively. 
The second term is the magnetic tension force (hoop-stress), which is always directed inward (in 
the — r direction). The third is the centrifugal force, which always acts in the +r direction, opposite 
to the hoop-stress. 

One key parameter of radial equilibrium in rotating jets is the radial index K of the magnetic 
pressure p^ = (5^ +5^)/2 gradient 

d log p* 
K=^-^^. (31) 

a log r 

Under the assumption that the magnetic field strength deceases outward near the jet surface, K 
should be less than zero. Furthermore, for specific cases of analytically asymptotic solutions for 
cylindrically symmetric axial flows with power-law radial distributions of the physical quantities, a 
lower limit on K must exist (Ostriker 1997). The magnetic field strength \B* \ [= (Bi+B^) ^'^] must 
decline with r more slowly than r^^ to ensure force balance and cylindrical coUimation. (That is, 
the magnetic field cannot have a gradient so steep that the corresponding magnetic pressure force 
will exceed the hoop-stress.) Taken together, these conditions require — 2 < K < 0. 

We now introduce another key parameter, the azimuthal (toroidal) Alfven Mach number: 

Ma,^^:^^, (32) 

where Va(|) = l^(|)|/v^- In the present paper, we have excluded the situation that jets are confined 
solely by the external hot ambient medium. So, to prevent radial expansion of the self-coUimated 
jets with K < 0, an upper-limit on M^i^ should be set: 

Ma0 < 1. (33) 

In the following sections, we will consider these two key parameters as diagnostics of radial 
force equilibrium in the rotating MHD jets. 



4.4.3. Non-rotating strongly magnetized jets 

If the jet has a low plasma-(3 (<^ 1) and is non-rotating (Ma(S^ = 0), the magnetic field in the 
jet will be in force-free equilibrium 
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That is, the gradient of the magnetic pressure will be balanced by the magnetic hoop-stress. Various 
types of the force-free equilibrium in jets have been investigated for different radial profiles of the 
magnetic twist $(r) (Appl, Lery, & Baty 2000; Lery, Baty, & Appl 2000). For astrophysical jets 
of magnetic origin, the azimuthal field is likely to be dominant, and the properties of the fastest 
growing CD kink (m = 1 ) mode become nearly independent of the details of the pitch profile (Appl, 
Lery, & Baty 2000). 



4.4.4. Rotating strongly magnetized jets 

If the jet is rotating {Mp^^ > 0), an equilibrium /orce-/ree magnetic field will not be possible. 
The inertia of the rotating plasma will play an important role in jet stability as M^i^ becomes large. 
This situation is quite different from the one that occurs in solar coronal loops. 

For an azimuthally sub-Alfvenic flow (0 < M^i^ < 1) in equilibrium, equation (30) will be 
satisfied. The hoop-stress is the only force acting inward and this must be balanced by the outward 
centrifugal and magnetic pressure gradient forces (k < 0). 

For an azimuthally trans-Alfvenic flow A/aij) ~ 1, the gradient of the magnetic pressure be- 
comes zero (i.e., K ~ 0), and force-equilibrium between the hoop-stress and the centrifugal force 
holds: 

-^ + ^^0. (35, 

r r 

If the radius r of the magnetized loop increases slightly, the centrifugal force will fall off faster 
than the hoop-stress, causing an inward-directed net magnetic tension that returns the loop to its 
original size. On the other hand, if r decreases slightly, the hoop-stress increase less rapidly than 
the centrifugal force, causing an outward-directed net centrifugal force that again returns the loop 
to its original size. Rapid rotation, therefore, has a potentially stabilizing effect on traditionally 
unstable twisted magnetic configurations. 



4.4.5. Numerical results for non-rotating and rotating jets: growth of axial modes 

We now inspect our numerical results more closely and analyze the stabilizing effects of 
rotation against the CD kink mode. We will continue to use the Fourier power spectrum analysis, 
applying it to the m = 1 mode and also extending it to growing modes in the axial (k) direction. We 
also will re-examine the distribution of 7^ and study a correlation between the maximum allowable 
magnetic twist 4>(r) and the toroidal Alfven Mach number Ma(\, for some selected parts of the jets. 
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Figure 17(1) shows a space-time (kj) diagramof the power spectrum of |/-''^(1,X)| for model 
A - 2. The shortest wavelength physically growing m=\ mode is Xmin ~ 1-0; perturbations with 
wavelength longer than this already have started growing hy t = 15.0. An enlarged snapshot of J^ 
for A-2at?= 19.5 (from Fig. 10) also is shown in Fig. 17 (2). Note the wiggled structure behind 
(upstream of) the R— S point with a destabilization length of ?i ~ 2.0. The 3rd panel of Fig. 17 (3) 
shows the transverse profiles from panel (2) of $(r) (for X = 2.0) and of Mai)) at two different axial 
positions in the jet (z=6.1 and 8.3). 

The value of $(r) at the point z = 6.1 increases gradually towards the inner jet. It exceeds the 
K-S critical value 4>crit = 271 for radial position \x\ < 0.4 and has a maximum (~ 10) near the jet 
axis. At z = 8.3, on the other hand, $(r) shows a sharp rise inward and has already exceeded ^cn\ 
by \x\ ~ 1.0. Between 0.4 < \x\ < 1.0, the magnetic twist is very large (> 20), followed by a sharp 
fall to around ~ 10. 

The values of Mp^^ at both axial points decrease gradually towards the inner jet, approaching 
close to within \x\ < 0.4. Because there is little rotation in the region \x\ < 0.4, the classical K-S 
criterion (25) can be used to analyze the stability of this shallow-atmosphere, mildly PFD case jet. 
Noting that the twist $(r) exceeds the K-S limit at z = 8.3 throughout most of the body of the jet, 
we conclude that the jet at that point is already unstable to the CD kink mode. On the other hand, 
because at z = 6.1 the K-S criterion is moderately violated in the region \x\ < 0.4, we conclude 
that the flow in this model is beginning to go unstable to the CD kink mode at that point. 

Figure 18 performs a similar analysis for another model from Fig. 10 — the shallow-atmosphere, 
highly PFD case A - 1. In this case, the shortest wavelength physically growing m = 1 mode 
changes with time, gradually shifting from X-aim ~ 2.5 to ~ 1.0. However, at the chosen model 
time {t = 21.5), perturbations with X < 2.0 have not yet started growing. [Compare Fig. 10 and 
Fig. 18 (1) and (2)]. Yet,inFig. 18 (3), the values of the magnetic twist $(r) at Z?o?/z axial positions 
(z = 6.1 and 8.3) appear to exceed the critical 4>crit for \x\ < 0.4 and have a maximum (~ 10) near 
the center. 

The reason for this seeming violation of the K-S stability criterion, and yet apparent jet sta- 
bility, can be seen in the M^i^ plot in Fig. 18 (3). The value of the azimuthal Mach number at 
both axial points is large for all radii in the jet, never dipping below Mai^ ~ 0.5. The large rota- 
tional velocity of the jet in this case provides a stabilizing influence over the classical K-S criterion 
throughout this jet region. 

Our nonlinear numerical results are consistent with the linear theory. A linear stability anal- 
ysis of the CD m = 1 mode was performed for a non-rotating cold (strongly magnetized) jet (in 
force-free equilibrium) by Appl (1996) and Appl, Lery, & Baty (2000). They showed that the 
most unstable wavenumber A:niax, as well as the maximum growth rate Im(a))niax, are only slightly 
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affected by the jet (axial) velocity Vj and the radial profile of the magnetic pitch rBz/B(^. (Note that 
the decrease of Im(a)) with increasing Ma = Vj/Va in their papers is associated with a decrease 
in Va, not an increase in Vy) In our results, unstable axial modes are excited over a large wave- 
length range. However, the wavelength X^nx corresponding to fc^ax seems to occur in the range 
2 < Xmax < 3 for both the super-Alfvenic (A - 2) and sub-Alfvenic (A - 1) cases [see Figs. 17 (1) 
and 18(1)]. Therefore, we do not see a strong correlation between kmax and Ma in our results, and 
this is qualitatively consistent with the linear analysis. 

These results lead us to conclude that the classical K-S criterion holds in approximately non- 
rotating jets (even the nonlinear regime), but that rotation in a jet creates a stabilizing effect against 
the CD kink mode. Rotating jets can remain stable even when the magnetic twist exceeds the 
classical K-S criterion. 



4.5. Sudden Destabilization of a Jet by Angular Momentum Loss 

Having shown that the V^ velocity component plays an important role in stabilizing the jet, 
we now investigate how the jet behaves if that rotation is suddenly removed. When (3 <^ 1, if Vc]) 
were to decrease gradually as the jet propagates, e.g., Visf0^z~^, then the magnetic field in the jet 
would shift quasi-statically toward a force-free (7 x B ~ 0) configuration. However, if part of the 
jet suddenly loses its toroidal velocity though some dynamical means, then that part will subject 
to a non force-free magnetic field (JxB^ 0). Here we investigate the sudden loss of rotation as a 
possible trigger process for the CD instability. 



4.5.1. The principal driving force of the CD kink 

First, we shall determine the principal driving force of the CD kink (m = 1) mode. Transverse 
profiles of the physical quantities in the intermediate stage at axial position z = 8.3 are plotted 
in Fig. 19. Panel (1) shows the specific centrifugal (Vl/r) and hoop-stress (— 55/r/p) forces for 
all four models. Panel (2) shows the distribution of magnetic pressure contributed by both the 
azimuthal and axial components p^ = (Bl+BJ)/!. 

In both models B - 1 and B - 2, the centrifugal force is balanced by the hoop-stress, so the 
flows are azimuthally trans-Alfvenic {Ma^ ~ 1). In addition, the magnetic pressure in these two 
models has quite a small gradient and so contributes little to the radial force balance. The flow is 
apparently CD stable at this axial position (z = 8.3), so we conclude that the azimuthal rotation has 
stabilized the flow. 
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In model A - 1, the centrifugal force is smaller than the hoop-stress, so the flow is azimuthally 
sub-Alfvenic (Ma^ < 1) and has a significant magnetic pressure gradient working with the centrifu- 
gal force to balance the hoop-stress. It is noteworthy, however, that the position of the magnetic 
pressure peak lies at the center of the jet (x, y)={0, 0): the force balance on either side of the jet is 
therefore symmetric. This situation is also apparently CD stable. 

In model A - 2, however, the relative contribution of the centrifugal force is essentially zero 
(Ma(\, — 0). The only forces trying to remain in balance are the magnetic hoop-stress and the radial 
magnetic pressure gradient. Figure 19, however, shows that these forces are not able to maintain 
a symmetric distribution as they try to balance each other. The hoop-stress at this axial position 
is asymmetric, and the (otherwise symmetric) magnetic pressure profile peak has been pushed 
to a position offset from the axis. The 3-D response of a force-free magnetic configuration to 
asymmetric hoop-stress, then, is not to create an asymmetric magnetic pressure profile but rather 
to push that symmetric profile off to one side. This asymmetry indicates that the CD modes have 
already begun to grow to a measurable level, and that this part of the jet is being accelerated in a 
negative radial direction {—x). 

We conclude that the CD instability is driven by asymmetric hoop-stress, created in a region 
of large magnetic twist, which cannot be balanced stably by a corresponding asymmetric magnetic 
pressure gradient. 



4.5.2. Origin of the very large magnetic twist 

We now examine the origin of this high magnetic twist in the jet. The twist is characterized 
by a relatively large \B^\ and a small B^. 

Figure 20 shows the intermediate-stage radial profiles (parallel to the r-axis along (|) = 0) of 
each component of the magnetic field {\B^\ and B^ and the magnetic pressure pj!^, as well as the 
initial distribution of the axial field component 5™' (same for all models). As we saw in Fig. 19 
(2), the distributions of the magnetic pressure for B - 1 and B - 2 are flat. This creates separate 
|5(|)| and B^ distributions that are low and high at different radii. 

For models A - 1 and A - 2 we see that, in the region 0.3 < r < 1.0, the magnetic pressure 
follows K ~ — 1 and K ~ — 2 distributions, respectively. Of particular interest in these plots is 
the behavior of B^ at these intermediate radii. Compared with the initial distribution, the axial 
field becomes stronger at small radii (r < 0.3) and remains about the same strength at large radii 
(r < 1.5). However, in the region 0.3 < r < 1.5 it dips substantially, particularly for model A - 2 
in Fig. 20. 
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We conclude that the very large magnetic twist (which leads to asymmetric hoop stress and 
a CD instability) is caused not so much by a large increase in B^ as by a large decrease in B^ at 
intermediate radii [see, especially, Fig. 17(3)]. 



4.5.3. Origin of the weak axial field causing the large twist: triggering of the CD kink instability 

by the reverse slow-mode shock 

Finally, we shall identify the physical mechanism that causes B^ to decrease and create the 
strong magnetic twist. In particular, this will occur when rotation is suddenly removed from the 
jet flow, triggering the pinch. We find that this occurs at a number of places in the jet flow, but is 
particularly strong right behind the R— S shock wave. 

Figure 21 shows in the intermediate evolutionary stage offset axial profiles (parallel to the 
z-axis, at (r, (|))=(0.5, 0)) of each component of the magnetic field and velocity. All models begin 
with the same power-law distribution of B^ with z [see panel (1)]. However, only model A - 2, 
and to a lesser extent A - 1, develop a large region of substantially low B^, especially in the region 
8.0<z<10. 

As pointed out in Sec. 4.2.2, there is a sudden rise in \B^\, V^, and V,^ around z = 10.0 for 
all models [see Fig. 21 (2), (3), and (4)], which corresponds to the passage of the F— F compres- 
sive/shock wave front. The models split into two classes [see panels (2) and (4)]: those with a 
large enhancement in \B(^\ and a small enhancement of Vc]) (A - 1 and A - 2) and those with a small 
enhancement in \B(^\ and a large increase in Vij, (B - 1 and B - 2). 

A large enhancement in \B^\, with a small enhancement in V,^, indicates that a considerable 
amount of energy has been converted from the toroidal (rotational) kinetic energy in the TAW to 
toroidal magnetic energy behind the F— F shock wave front. 

The sudden drop of both \B^\ and V,^atz = 9.2 in model A - 2 (dashed lines) and at z = 5.5 for 
B - 2 (dotted lines) occurs with the passage of the F— S shock wave front. We also can identify the 
R— S shock wave front near axial position z = 8.5 for model A - 2. At this point \B^\ is enhanced 
once again (panel (2)), but V(^ has a sharp dip before rising to its injected value for z < 8.0. The 
contact discontinuity lies in between the two slow-mode shock waves — at z ~ 8.7. 

The behavior of these quantities at the R— S shock wave is especially interesting. While \B^\ 
simply increases, the dip in V,^ has important consequences. In this transition, part of the jet near 
z = 8.3 loses toroidal velocity (angular momentum) but still gains enhanced toroidal magnetic field. 
The centrifugal support against the hoop-stress, therefore, is lost, causing a sudden squeezing of 
of the poloidal magnetic flux \Bp\ = [B^ + B^y^^ towards the central (z) axis. As a result, a 
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concentration of 5^ flux occurs at the jet center (r < 0.3), at the expense of B^ flux at intermediate 
radii (0.3 < r < 1.5) (see also A - 2 in Fig. 20). 

To summarize, we believe that the CD kink mode is triggered as follows. The process begins 
when rotating, magnetized jet material experiences a sudden loss of kinetic angular momentum to 
the magnetic field. This destroys the (stable) balance between centrifugal force and hoop-stress, 
producing a strong pinch of B^ flux toward the z axis, depleting B^ flux in the intermediate radial 
region of the jet, and further enhancing the magnetic twist. With no rotation to stabilize this large 
twist, the K-S criterion is suddenly strongly violated and the jet becomes unstable to the CD kink 
mode. Any slight radial asymmetry in the hoop-stress produces a helical kink in the current, with 
no stabilizing restoring force. This process can occur at various places in the jet flow, but appears 
strongest right behind the R— S shock wave, where there can be a pronounced dip in the rotational 
velocity. We believe that this is one of the most promising physical processes for producing a CD 
kink instability in rotating jets. 



4.6. Suppression of MHD KH Instabilities due to the External Magnetized Winds 

We now show how the helically magnetized winds that surround the core of the jet can sup- 
press MHD KH instabilities. 

The modified stability criterion for MHD KH instabilities is as follows (Hardee & Rosen 
2002): 

Ay2_y2^<0, (36) 

where AV = Vj — Ve is the velocity shear between the jet and the external medium, and Vas = 
[(pj + pe)(5?+5g)/(pjPe)]^'^ is the surface Alfven speed. Inequality (36) shows that a magnetic 
field and an outgoing wind in the external medium act to stabilize the flow against MHD KH 
surface modes. Figure 22 shows, for all four of our simulations, the transverse profiles in the 
X direction of the difference between AV^ and V^^, as well as the total Alfven Mach Number 
Ma (= IVI/Va)- Two different evolutionary stages are compared: the intermediate stage at axial 
position z = 8.0 and the final stage at z = 14.0. 

In the intermediate stage, the MHD flow becomes trans- to super-Alfvenic for model A - 2, 
but remains sub-Alfvenic for all other models, as shown in Fig. 22 (1-2). Moreover, for model 
A - 2 Ma varies smoothly at a radius (r ~ \x\ ~ 0.4), roughly the boundary between the current- 
carrying jet and the outer magnetized wind (see also Fig. 1 1). As is seen in Fig. 22 (1-1), inequality 
(36) is satisfied for all models, independent of the absolute value of Ma. (Note that the distribution 
for A - 2 is already slightly asymmetric, even though the flow is stable to MHD KH instabilities. 
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This indicates that the helical kink that develops must be caused by not the MHD KH instability 
but another possible mechanism — the MHD CD instability in this case.) 

In the final stages, the flow of jet still remains sub-Alfvenic for model A - 1, but trans- to 
super- Alfvenic for all other models, as shown Fig. 22 (2-2). By this time all of these current- 
carrying jet systems, except model B - 2, have been distorted by the CD instability. As a result, the 
distributions of Ma for the three distorted models are decidedly asymmetric. Note that, whether the 
instability has grown (A - 1, A - 2, and B - 1) or not (B - 2), inequality (36) still holds everywhere 
for these models. 

We conclude that the MHD KH surface modes are effectively suppressed, because the external 
magnetized winds reduce the velocity shear in the transverse direction. Our MHD jets, therefore, 
are distorted not by KH instabilities but pure CD instabilities, even when the flow itself becomes 
super- Alfvenic. 



4.7. Entrainment of Ambient Material in PFD Jets 

In Fig. 23 we plot the axial profiles of the density and the ratio of Alfven and sound speeds 
in the final stage (similar to Fig. 1). Note that the distributions of mass entrainment and Va/Cs(~ 
P 1/2) depend not on the initial conditions (initial gradients in the ambient medium), but on the 
nature of the flow itself (on whether it is strongly Poynting flux dominated or not). As the jets 
propagate through the ambient medium, they adjust themselves to the decreasing external pressure. 
If the ambient medium has no significant field strength ((3 ^ 1), the jets expand, the density ratio 
Pj/Pe becomes smaller, and the jets are rarefied. In this case, the gradient of the external medium 
may play a crucial role in the entrainment of matter. 

However, in the opposite case (|3 ^ 1), the jets have considerable field strength and are pre- 
vented from expanding by the self-confining hoop-stress. This means the external density gradient 
might not have as much influence on mass entrainment. 

The mass entrained in the mildly PFD flows is much greater than in the highly PFD flows. In 
particular, p decreases slower than ~ z^^ for both models A - 2 and B - 2, but roughly as ~ z^^ for 
models A - 1 and B - 1 as seen in the Upper panel of Fig. 23. This is partly due to the difference 
in injected mass flux into the computational domain. (Low mass flux implies high Poynting flux.) 
However, some of the effect also appears to be related to the initial ambient density gradient, with 
a steep atmosphere causing a higher mass entrainment in jets than the slowly decreasing case. The 
highest entrainment occurs for model B - 2, where the flow is both mildly PFD and the atmosphere 
has an initially steep density gradient. An enhanced mass distribution also can be seen in model B 
- 1, even though it is highly Poynting flux dominated. On the other hand, for model A - 1, which 
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neither is mildly PFD nor has an initially steep density gradient, essentially no enhancement (in 
fact, a slight decrease) can be seen. Denser jets in general are predicted to remain stable against 
MHD KH instabilities beyond the Alfven surface and have been found to be more robust than their 
less dense counterparts (Hardee, Clarke, & Rosen 1997). 

The Va/Cs distributions also depend strongly on the nature of the flow as seen in the Lower 
panel of Fig. 23. Compared with their initial distributions of Va/Q^ all models except A - 2 
become significantly flattened in the z direction. Only this shallow-atmosphere, mildly PFD model, 
develops an even steeper distribution |3^^/^. 

Even with this mass entrainment, however, our jets still keep their highly magnetized state 
^a/Cs ^3 (P < 0.1) during their dynamical evolution, except in the region right behind the F— F 
compressive/shock wave front where compression is high. 



4.8. Saturation and Advection of Kinked Jet Structures 

As we have seen in Fig. 15 and 16, the patterns created by the CD instability can be advected 
by the propagating jet even after saturation of that instability has occurred. This means that the CD 
instability does not destroy the interior of the jet. Rather, it only distorts the interior into a semi- 
permanent wiggled structure. In the final stage of the dynamical evolution, the jets still remain 
strongly magnetically dominated (|3 < 0.1), and the flow never attains a state of MHD turbulence 
(which requires |3 ^ 1). If the magnetic diffusion time scale {i.e., dissipation time scale of the 
current) is much longer than the jet propagation time scale for trans-Alfvenic flow (the Alfven 
transit time scale) (Koide, Nishikawa, & Mutel 1996), these patterns will persist for some time as 
the flow advances. And the plasma flow itself will appear to travel in a true helical pattern as it 
follows the magnetic backbone of the helix (see the velocity vectors in Figs. 12 and 16). 



5. CONCLUSIONS AND APPLICATION 

5.1. Summary and Conclusions 

By performing 3-D non-relativistic MHD simulations we have investigated in detail the non- 
linear dynamics of PFD jets and the development of asymmetric instabilities. The PFD jets, pow- 
ered by TAWs, propagate into extended stratified atmospheres. The presence of various compres- 
sive MHD waves in jets, which can steepen into MHD shock waves, plays an important role in 
the nonlinear evolution of the system. Effects such as shock waves or re-distribution of the ax- 
ial current profile in the post-shocked region strongly affect the excitation of the CD instabilities. 
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The growth of the CD m = 1 mode causes a kinked, 3-D spatial helical structure, which may be 
responsible for the wiggled structures seen in images of some observed jets. 

We have performed four simulations with different combinations of the following initial pa- 
rameters: highly (Fexb/Fioi ~ 0.9) or mildly (Fexb/Fioi ~ 0.6) PFD jets with decreasing or con- 
stant Va external ambient medium along the central axis (z). From these simulations we conclude 
the following: 

1 . Due to a centrifugal effect, rotating jets can be stabilized against the CD kink mode beyond 
the point predicted by the classical K-S criterion. Radial force balance in strongly magne- 
tized (C^ <^ V^) rotating jets, therefore, can involve a centrifugal force component; it need 
not involve only a pure electromagnetic force-free field. 

2. Non-rotating jets will be subject to the CD kink mode when the classical K-S criterion is 
violated. The kink (m = 1) mode grows faster than the higher order modes (m > 1). 

3. The driving force of the CD kink mode is an asymmetrically distribution of hoop-stress 
(the magnetic tension force). This is caused by a sudden decrease of jet rotation and a 
concentration the poloidal magnetic flux toward the central axis via nonlinear processes that 
occur at a reverse slow-mode MHD shock wave. As a result, the magnetic twist become 
larger than the critical value specified by the K-S criterion. 

4. The linear growth rate of the CD kink mode does not depend on the flow speed. Instead, it 
grows on time scales of order of the local Alfven crossing time. Kinetic-flux dominated jets 
that have a large Alfven Mach number Ma (~ Mfm in C^ <^ V^) are not subject to the CD 
kink mode, because the jet propagation time is shorter than the CD instability growth time. 

5. We do not see growing surface modes of the MHD KH instability even when our flows 
become super- Alfvenic. The naturally-occurring, external helically magnetized wind, which 
is (quasi-) axially current-free (7^ ~ 0), surrounds the well-coUimated current-carrying jet. 
This wind reduces velocity shear between the jet and the external medium, suppressing the 
growth of KH surface modes. Therefore, we have been able to investigate exclusively the 
nonlinear behavior of pure CD instabilities. 

6. Under strongly magnetized ambient conditions, the amount of mass entrained and the ratio 
Va/Q in the jet depend critically on the nature itself of the flow (/. e. , on the level of Poynting 
flux domination) and much less on the external density or gas pressure gradients. The prop- 
erties of jets observed thousands of Schwarzschild radii from the galactic center, therefore, 
may nevertheless carry important clues about the nature of the central engine itself. 
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The results presented in this paper are valid in the non-relativistic regime. For relativistic 
motion {V ~ c), the electric field E may be comparable to the magnetic field B, but it will never 
dominate (B^ — £^ > in the force-free limit). As in these non-relativistic simulations, therefore, 
the CD instabilities will not propagate in the rest frame of the jet (jet co-moving frame). If we 
choose the inertial frame to be the rest frame of the jet fluid, E will vanish for a perfectly conducting 
fluid (the MHD condition). With E (and the displacement current D) negligible, the rotation of a 
relativistic PFD jet (in which the energy density of the electromagnetic field is much greater than 
that of the particles) still will play an important role in stabilizing the flow against the growing CD 
instability, even in the relativistic limit. Full treatment of nonlinear CD instabilities in relativistic 
PFD jets will be given in a forthcoming paper. 



5.2. The CD Instability as the Origin of Kinked Structure in AGN and Pulsar Jets 

The behavior of our simulations is consistent with helical flows seen in such sources as 3C 
345 (Zensus et al. 1995) and 3C 120 (Gomez et al. 2001). Helical structures appear to persist 
for longer than component propagation times, and the plasma components themselves appear to 
propagate on helical trajectories. 

Of course, a three-dimensional helical shape in a real jet could be produced by precessional 
motion at the jet origin instead of by asymmetric CD (m > 0) modes. In both cases the center of 
the actual jet flow would deviate from the jet central axis during dynamical evolution. However, 
precessional jets do not show true helical motion as their components move outward from the 
galactic center. They simply move ballistically, with the turning radius becoming larger as the jet 
moves further from the origin. 

In the case of asymmetric KH surface modes, such a deviation of the flow from the jet axis 
is not even possible (see, e.g.. Fig. 5 in Hardee & Stone 1997). While it might be so for the 
asymmetric KH "body" modes, the growth range for these modes is limited (see, 2.4.1). 

Therefore, if most parsec-scale helical jets in AGN are confirmed to be true helical flows, 
then the CD instability in a Poynting flux dominated jet may be a viable model for producing the 
kinked, helical structures seen in these objects. An excited wave disturbance, such as occurs in 
a KH mode, may not always be needed in order to explain the kinked morphological features in 
these jets. 

High angular resolution, time-dependent X-ray observations by the Chandra X-ray observa- 
tory displayed year-scale variations in jet morphology in the Vela (Pavlov et al. 2003) and Crab 
(Mori et al. 2004) pulsar wind nebulae. These sources showed a growth of their overall kinked 
structure as the jet moves outward. The Crab jet is 10 times larger and varies 10 times more slowly 
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than the Vela jet. (Time scales and jet widths are 150-500 days and 2.9 x 10^^ cm for the Crab, and 
10-30 days and 3.0 x 10^^ cm for Vela, respectively.) These kinked patterns seem to be co-moving 
with the main flow, and the time scales reported are proportional to the Alfven crossing time (with 
Va being approximately the speed of light for these ultra-relativistic plasmas in equipartition with 
the magnetic field). All of these features are consistent with our results on CD instabilities in PFD 
coUimated flows. Pulsar jets, therefore, are one of the strongest cases for a Poynting flux dominated 
astrophysical flow that has been distorted by current-driven instabilities. 
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in part, by a NASA Astrophysics Theory Program grant. This research was performed at the Jet 
Propulsion Laboratory, California Institute of Technology, under contract to NASA. Numerical 
computations were carried out on Fujitsu VPP5000/32R at the Astronomical Data Analysis Center 
of the National Astronomical Observatory, Japan. Visualizations were performed on SGI Origin 
2000 at the Supercomputing and Visualization Facility of Jet Propulsion Laboratory and on other 
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A. THE INITIAL FIELD CONFIGURATION 

The initial current-free field is explicitly given by the following formula: 

B= (Bx.By.B^) = (5^cos(|)-5(|,sin(|), 5;-sin(|)4-5(|)Cos(|), 5j;), (Al) 
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On the z-axis (r = 0) we have 
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Parameters r„ and z„ are the radius and position of the nth loop (ring) current, respectively, which 
are set as follows: r„ = 1.0 for all n, and the positions of pair of current loops are at Zn = —3.0 and 

2Zmax + 3.0. 
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Fig. 1 . — Initial profiles along the z-axis of (1) the axial component of the magnetic field B^ for all 
models and (2) density p and ratio of Alfven velocity to sound velocity Va/Q in logarithmic scale 
for models A (light gray solid line) and B (dark gray solid line). The broken line indicates a z^^ 
power law variation for comparison. 
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Fig. 2. — Time variation of the energy fluxes injected though the z = plane. Kinetic energy 
flux F]jin and Poynting flux Fexb are normalized by the total energy flux Ftot = ii^in +F£xB +^th 
(thermal energy flux). Both Case 1 {Highly PFD Jets) and Case 2 {Mildly PFD Jets) are shown. 
The solid line shows Model A {slowly decreasing ambient density) and the broken line Model B 
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Fig. 3. — Time evolution of the region (\x\ < 2.25, 0.0 < z < 18.0) for Model A-1. Contour 
images of the density distribution (logarithmic scale) are shown along with the poloidal velocity 
field in the x — z plane. Model times are given at upper left in each of the four panels. The length 
of the arrow at upper right in each panel shows the unit velocity. 
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Fig. 4. — Similar to Fig. 3, but for Model A - 2. 
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Fig. 5. — Similar to Fig. 3, but for Model B - 1. 
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Fig. 6. — Similar to Fig. 3, but for Model B - 2. 
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Fig. 7. — Offset axial velocity profiles parallel the z-axis at {x, y) = (—0.025, 0.0) in the interme- 
diate stage for each of the models. Shown are the components of velocity (V^, Vy, V^) and the local 
Va and Q at model times ? = 21.5 (Upper-Left: A- l),t = 19.5 (Upper-Right: A-2),t = 10.0 
(Lower-Left: B - 1), and t = 10.0 (Lower-Right: B - 2). The labels F— F, F— S, and R— S indicate 
the positions of the forward fast-mode, forward slow-mode , and reverse slow-mode compressive 
MHD (magnetosonic) waves, respectively. Some of these will steepen into MHD shock waves dur- 
ing the dynamical evolution. Each time selected corresponds to when the F— F compressive/shock 
wave fronts reach approximately z = 10.0. 



53 



•A-1: t=£l.5 



H' « 



I 



-1 f — 



H^ 



1 



^^^; J 



f 



h F-F 



-1 1 1 r 



— 1 — 



14 



te sS 



h o 



1^ o i- 



W^ 



iS. iG 




Fig. 8. — Snapshots of the azimuthal velocity Vy in the intermediate stage for each of the models. 
The labels F— F, F— S, and R— S are the same as those in figure 7. 
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Fig. 9. — Similar to figure 8, but for By. 
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Fig. 10. — Snapshots of the axial current density J^ in the intermediate stage for each of the models. 
The figures show a closed circulating current system, in which one current path occurs near the 
central (z) axis and co-moves with the PFD jet (J^^), while the other conically-shaped current 
returns outside (J^^). 
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Fig. 1 1 . — Transverse profiles parallel to the ;c-axis of the force-free parameter Xff (Upper) and of 
the negative of the axial current density — 7^ (Lower) in the intermediate stage. Two axial positions 
are shown: z = 2.2 (Left) and z = 8.3 (Right). |?iff | = 1 indicates that J is parallel or anti-parallel to 
the local B direction (i.e., force-free). The reversal of the sign of Xff from -1 to 1 (seen in the Upper 
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seen in both Lower panels). 
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Fig. 12. — Snapshots of the axial current density J^ in the final stage for each of the models. 
Each time selected corresponds to when the F— F (A-l and A - 2) or the F— S (B-1 and B 
- 2) shock wave fronts approach or exceed approximately z — 18.0. Note that these figures are 
2-D meridional slices only through the 3-D jet structure. While the distributions of /J'^ appear 
asymmetric in models A - 1, A - 2, and B-1, the actual jet structures in these models are a 
three-dimensional spiral helix. 
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Fig. 13. — Similar to figure 7, but for the final stage in each of the models. 
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Fig. 14. — Growth of unstable CD modes for each model sequence, showing the time variation of 
the azimuthal Fourier power in modes m = 1 — 4, on a natural log scale. Solid lines are derived by 
fitting the slopes of the m = 1 mode growth to a linear function. 
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Fig. 15. — Space-time diagram showing the growth of the ratio of transverse and axial magnetic 
field —By/B^ during the interval t = 0.0 — 45.0. Values are measured near, and parallel to, the z 
axis (x = -0.025, y = 0.0) for Model A - 1. 
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Fig. 16. — Snapshots of the radial specific power generated by the Lorentz force (7 x fi)^y,./p in 
the final stage for each of the models. The r-component of the Lorentz force is composed of both 
(i) the magnetic pressure gradient of the azimuthal and axial components —{d/dr)[{Bi +Bf)/2] 
and of (ii) the hoop-stress —Bl/r. Positive power {Red) indicates the net acceleration (increase 
in radial velocity); negative power {Blue) indicates net deceleration. There is a strong correlation 
between the power distribution and the wiggled structures seen in models A - 1, A - 2, and B - 1. 
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Fig. 17. — Analysis of the K-S criterion for Model A - 2 in the intermediate stage. (1) Space-time 
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Fig. 18. — Similar to Fig. 17, but for Model A - 1. 
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field components B^, \B^\, and the magnetic pressure p^ (logarithmic scale) in the intermediate 
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Fig. 21. — Offset axial profiles, parallel to the z-axis, at (r, (|)) = (0.5, 0) in the intermediate stage 
for each of the models. (1) Axial magnetic field component B^ (the initial distribution is also 
plotted by light gray thick solid line). (2) Absolute value of azimuthal magnetic field component 
|5(|)|. (3) Axial velocity V^. (4) Azimuthal velocity Vh,. 
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Fig. 22. — Transverse profiles in the x direction of velocity shear and Alfven Mach number in the 
intermediate stage (Left panels at axial position z = 8.0) and final stage (Right panels at z = 14.0) 
of evolution. Upper panels: the difference between the square of the velocity shear AV^ and the 
square of the surface Alfven velocity V^^. Lower panels: the total Alfven Mach number Ma. The 
shaded area in the Upper figures shows the MHD KH unstable region (AV^ > V^s)- ^'^ the first 
panel, models B - 1 and B - 2 do not appear, because their velocity difference lies below —1. 
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Fig. 23. — Axial profiles of (1) the density p {Upper) and (2) the ratio of Alfven velocity to sound 
speed Va/Q (Lower) in the final stage of evolution. Each thick solid line represents the initial 
(t = 0.0) distribution for Model A (light gray) and Model B (dark gray). 



69 



Table 1. Units of Physical Quantities for Normalization. 



Physical Quantity Description Normalization Unit 

t .^^^_^^^^^_. Time Xao (= ^o/^Ao) 

L (= ^Jx^+y^ + z^) Length Lq 

p Density po 

p Pressure Po^AO 

V Velocity Vao 



B Magnetic field •, AnpaVjj^ 

J Current Density J AnpoVJj^/Lo 



Note. — The initial value of the density po and Alfven veloc- 
ity Vao at the origin (x, y, z) = (0, 0, 0) are chosen to be the typical 
density and velocity in the system. That is, the initial dimension- 
less density p' and Alfven velocity V^ at the origin are set to unity. 
The unit length scale Lo is approximately the same as the initial jet 
diameter at z = 0. So, we can define a characteristic time scale, the 
initial Alfven crossing time Tao at the origin, which is associated 
with Lo and Vao- The time is normalized with Tao> so the dimen- 
sionless Alfven crossing time also is set to unity: T^ = 1. 



Table 2. Model Parameters and Flow Properties. 



Model a" r'' VQ'' Wa'^ Vp'= Fexb/F,„ 



'0.6 
'0.9 
'0.6 
'0.9 



"Power index of p o= 1^1". 

''Polytropic index of p « p"^. 

'^Gradient of the sound velocity along the z-axis. 

''Gradient of the Alfven velocity along the z-axis. 

■^Gradient of the plasma-P along the z-axis. 

'Time averaged Poynting flux FexB injected into the system 
though the z = plane, normalized by the total energy flux Ftoi- 

Note. — Models A represent slowly decreasing atmospheres 
with decreasing Va. Models B represent steeply decreasing at- 
mospheres with constant Va- The signs "\", "/"', and "^" 
indicate negative, positive, and no gradients along the z-axis, re- 
spectively. 
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Table 3. Comparison of Growth Rates for All Models. 



Model 


Im((n) 


Growing 


azimuthal modes (m) 


Saturation 


A-1 .. 


0.31 




1,2 


No 


A-2 .. 


0.54 




1 


Yes 


B-1... 


0.51 




1 


No 


B-2... 


0.42 




1 


Yes 



